Introduction to Transparent Bayesian Data Analysis (chi22-course)
Abhraneel Sarma and Fumeng Yang
08 March, 2022
- Introduction
- Dataset 1
- Model 1. Wilcoxon signed rank test
- Model 2. Bayesian t-test
- Step 1: Build the model (likelihood)
- Step 1: Build the model (priors)
- Step 1: Build the model (prior predictive checks)
- Step 2: Computing the posterior probability
- Step 3: Evaluate the fit
- Step 3: Evaluate posterior predictions
- Step 3: Evaluate posterior predictions
- Step 3: Compare posterior predictions
- Compared to the frequentist estimates
- Model 3. Ordinal Logistic Regression
- Dataset 2
Introduction
In this document, we will outline the Bayesian analogs of the statistical analyses described here (and in the previous course lecture).
Helper functions.
# this function extracts results from different models and generate results of the same format to be used in visualizations
tidy.wrapper = function(model) {
if (class(model) == "lm") {
tidy(model, conf.int = TRUE) %>%
select(-c(statistic, p.value)) %>%
mutate(model = "Frequentist") %>%
select(model, everything())
} else if (class(model) == "brmsfit") {
tidy(model) %>%
filter(effect == "fixed") %>%
select(c(term:conf.high)) %>%
mutate(model = "Bayesian") %>%
select(model, everything())
} else {
stop("unknown model class")
}
}Dataset 1
Let’s first load and look at the data in a table
dataset = readr::read_csv("02-bayesian_stats/data/blinded.csv")
head(dataset)| experiment | condition | effectiveness |
|---|---|---|
| 1 | no_graph | 7 |
| 1 | graph | 6 |
| 1 | no_graph | 9 |
| 1 | no_graph | 4 |
| 1 | graph | 6 |
| 1 | no_graph | 7 |
dataset %>%
mutate(effectiveness = fct_rev(factor(effectiveness, levels = 1:9)),
experiment = as.factor(experiment)) %>%
# stacked bar plot
ggplot(aes(x = condition, fill = effectiveness)) +
geom_bar(position = "stack", stat="count") +
# plot data for different experiments as small multiples
facet_wrap(. ~ experiment) +
# grey color scale is robust for colorblind
scale_fill_brewer(palette="Purples", drop = FALSE) +
# horizontal plot
coord_flip() +
# legend
guides(fill = guide_legend(reverse = TRUE)) 
As we can see above, the original dataset contains results from four different experiments. For the purposes of this lecture, we will confine ourselves to the first experiment.
exp1.data = dataset %>%
filter(experiment == 1)
head(exp1.data)| experiment | condition | effectiveness |
|---|---|---|
| 1 | no_graph | 7 |
| 1 | graph | 6 |
| 1 | no_graph | 9 |
| 1 | no_graph | 4 |
| 1 | graph | 6 |
| 1 | no_graph | 7 |
Model 1. Wilcoxon signed rank test
This is a non-parametric test which we will skip for now. Although, there exists Bayesian non-parametric methods, they are more advanced for this lecture.
Model 2. Bayesian t-test
This is the linear model equivalent for the paired sample t-test.
dataset1.lm.freqt <-
lm(
effectiveness ~ condition - 1, # we remove intercept, which is equivalent
data = exp1.data
)Step 1: Build the model (likelihood)
Although we know that we are going to use a t distribution, we still plot the data to set how it looks like. Usually, this will help us decide our likelihood function.
exp1.data %>%
ggplot(., aes(x = effectiveness)) +
geom_histogram(fill = '#351c75',
color = NA,
binwidth = .5)
Then we write the formula using brm formula (bf)
function.
dataset1.brm.bayesiant.formula <- bf(effectiveness ~ condition - 1,
family = student())Step 1: Build the model (priors)
We check the priors that we should set
The blank ones are flat (uniform) priors. These are improper priors
and usually needed to be set. We can and should adjust other priors
given by brm.
as_tibble(get_prior(dataset1.brm.bayesiant.formula, data = exp1.data))| prior | class | coef | group | resp | dpar | nlpar | bound | source |
|---|---|---|---|---|---|---|---|---|
| b | default | |||||||
| b | conditiongraph | default | ||||||
| b | conditionno_graph | default | ||||||
| gamma(2, 0.1) | nu | default | ||||||
| student_t(3, 0, 2.5) | sigma | default |
Our priors.
dataset1.brm.bayesiant.priors = c(
# there's a lot of data so even fairly "strong" priors are going to not matter so much here
prior(normal(0, 1), class = "b"),
prior(student_t(3, 0, 1), class = "sigma"))Step 1: Build the model (prior predictive checks)
We do a prior predictive checks to see the predicted range if we only use priors. This helps us understand the predicted range of priors and formula, and also helps us understand if our priors and formula are appropriate.
Actually, the predicted range is quite wide for this model. We could consider shrink priors.
dataset1.brm.bayesiant.priorchecks <-
brm(
dataset1.brm.bayesiant.formula,
prior = dataset1.brm.bayesiant.priors,
data = exp1.data,
backend = "cmdstanr",
sample_prior = "only",
file = "02-bayesian_stats/rds/dataset1.brm.bayesiant.priorchecks.rds"
)
dataset1.bayesiant.yprior <-
posterior_predict(dataset1.brm.bayesiant.priorchecks)
ggplot() +
geom_density(
aes(x = dataset1.bayesiant.yprior),
color = '#351c75',
alpha = .5,
size = 1
) +
xlab('prior draws')
Step 2: Computing the posterior probability
dataset1.brm.bayesiant =
brm(
dataset1.brm.bayesiant.formula,
prior = dataset1.brm.bayesiant.priors,
data = exp1.data,
# You may not need this if rstan works with your computer.
backend = "cmdstanr",
# There are many other parameters you can play with.
# Here we use default parameters to simplify the lecture.
# We save the model
# But if you have this file, brm only loads the model for this file
file = "02-bayesian_stats/rds/dataset1.brm.bayesiant.rds"
)Step 3: Evaluate the fit
summary(dataset1.brm.bayesiant)## Family: student
## Links: mu = identity; sigma = identity; nu = identity
## Formula: effectiveness ~ condition
## Data: exp1.data (Number of observations: 123)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 6.63 0.18 6.28 6.97 1.00 3934 2469
## conditionno_graph 0.15 0.24 -0.31 0.63 1.00 3929 2920
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.29 0.12 1.05 1.52 1.00 2543 2387
## nu 16.63 11.78 4.16 46.91 1.00 2673 2496
##
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Step 3: Evaluate posterior predictions
We then call bayesplot to draw posterior distributions
and MCMC traces for each parameter.
We should see that the four chairs are mixed well, meaning that the model converged.
color_scheme_set(scheme = "purple")
plot(dataset1.brm.bayesiant)
Step 3: Evaluate posterior predictions
We draw from posterior predictive distributions. They are many ways
you can do this. Here we use posterior_predict from
brms.
dataset1.bayesiant.y <- exp1.data$effectiveness
dataset1.bayesiant.yrep <-
posterior_predict(dataset1.brm.bayesiant)
head(as_tibble(dataset1.bayesiant.yrep[, 1:11])) | V1 | V2 | V3 | V4 | V5 | V6 | V7 | V8 | V9 | V10 | V11 |
|---|---|---|---|---|---|---|---|---|---|---|
| 5.874652 | 5.619724 | 7.546285 | 3.835248 | 6.821533 | 7.522359 | 8.579928 | 8.501289 | 5.921710 | 5.138925 | 4.043591 |
| 9.352873 | 5.289061 | 7.801147 | 7.808275 | 8.707273 | 5.872816 | 6.769362 | 6.391579 | 7.907620 | 8.034020 | 9.689779 |
| 8.230550 | 6.520340 | 9.969986 | 6.113050 | 8.303660 | 5.467097 | 5.144298 | 5.318557 | 7.244001 | 5.623543 | 6.600704 |
| 5.990564 | 7.715373 | 4.356735 | 6.991714 | 9.004902 | 8.646328 | 5.914456 | 1.310394 | 8.073306 | 6.469923 | 6.221491 |
| 10.562348 | 4.278225 | 8.714006 | 6.856854 | 8.130915 | 8.076192 | 6.941659 | 5.471274 | 5.399386 | 5.139248 | 8.347167 |
| 6.748229 | 7.159318 | 8.387614 | 6.027675 | 4.142546 | 7.984125 | 8.204699 | 6.176600 | 7.854791 | 4.639312 | 5.572873 |
We use bayesplot to plot the first eight draws of posterior predictions and the original data. Again, there are many ways to do this (more in Lecture 3).
ppc_hist(y = dataset1.bayesiant.y,
yrep = dataset1.bayesiant.yrep[1:8,],
binwidth = .5)
ppc_dens_overlay(y = dataset1.bayesiant.y,
yrep = dataset1.bayesiant.yrep[1:30,])
ppc_stat_grouped(y = dataset1.bayesiant.y,
yrep = dataset1.bayesiant.yrep,
group = exp1.data$condition)
Step 3: Compare posterior predictions
We compare the two conditions using conditional distributions. add_epred_draws from tidybayes is the most convenient way…
dataset1.bayesiant.posterior_epred <-
tibble(condition = c('graph', 'no_graph')) %>%
add_epred_draws(dataset1.brm.bayesiant,
re_formula = NA,
allow_new_levels = FALSE) %>%
ungroup()
kable(head(dataset1.bayesiant.posterior_epred))| condition | .row | .chain | .iteration | .draw | .epred |
|---|---|---|---|---|---|
| graph | 1 | NA | NA | 1 | 6.55796 |
| graph | 1 | NA | NA | 2 | 6.77053 |
| graph | 1 | NA | NA | 3 | 6.61823 |
| graph | 1 | NA | NA | 4 | 6.61574 |
| graph | 1 | NA | NA | 5 | 6.35145 |
| graph | 1 | NA | NA | 6 | 6.46093 |
We then transform the data, subtract the two conditions from each other, and calculate the credible intervals.
dataset1.bayesiant.posterior_comparison <-
dataset1.bayesiant.posterior_epred %>%
select(-c(.chain, .iteration, .row)) %>%
compare_levels(variable = .epred, by = condition)
# just use head(dataset1.bayesiant.posterior_comparison) if you don't need to knit n html
head(dataset1.bayesiant.posterior_comparison)| .draw | condition | .epred |
|---|---|---|
| 1 | no_graph - graph | 0.1566460 |
| 2 | no_graph - graph | 0.1408930 |
| 3 | no_graph - graph | 0.0747382 |
| 4 | no_graph - graph | 0.3541770 |
| 5 | no_graph - graph | 0.7305980 |
| 6 | no_graph - graph | 0.3816200 |
dataset1.bayesiant.posterior_comparison %>%
mean_qi(.epred)| condition | .epred | .lower | .upper | .width | .point | .interval |
|---|---|---|---|---|---|---|
| no_graph - graph | 0.1544225 | -0.3123828 | 0.627965 | 0.95 | mean | qi |
We can also plot the credible interval. Lecture 3 will cover more.
dataset1.bayesiant.posterior_comparison %>%
median_qi(.epred) %>%
ggplot() +
geom_point(aes(x = .epred, y = condition), color = DARK_PURPLE, size = 3) +
geom_errorbarh(
aes(xmin = .lower, xmax = .upper, y = condition),
color = DARK_PURPLE,
alpha = .5,
size = 2,
height = 0
) +
geom_vline(aes(xintercept = 0), linetype = "dashed", color = "gray") +
coord_cartesian(ylim = c(0, 2), xlim = c(-1, 1)) +
xlab('') + ylab('')
Compared to the frequentist estimates
bind_rows(tidy.wrapper(dataset1.lm.freqt),
tidy.wrapper(dataset1.brm.bayesiant)) %>%
ggplot() +
geom_pointrange(
aes(
x = model,
y = estimate,
ymin = conf.low,
ymax = conf.high,
color = term
),
position = position_dodge(width = 0.2)
) +
scale_color_brewer(palette = "Set1") +
ylab('effectiveness') +
scale_y_continuous(breaks = 1:9, limits = c(1, 9)) +
coord_flip()
Model 3. Ordinal Logistic Regression
A more appropriate to analyze ordianl data is to use ordinal logistic regression. In this section, we reanalyze the data using Bayesian rdinal logistic regression.
Step 1: Build the model (likelihood)
dataset1.brm.olr.formula <-
bf(effectiveness ~ condition,
family = cumulative("logit"))Step 1: Build the model (priors)
as_tibble(get_prior(dataset1.brm.olr.formula, data = exp1.data))| prior | class | coef | group | resp | dpar | nlpar | bound | source |
|---|---|---|---|---|---|---|---|---|
| b | default | |||||||
| b | conditionno_graph | default | ||||||
| student_t(3, 0, 2.5) | Intercept | default | ||||||
| Intercept | 1 | default | ||||||
| Intercept | 2 | default | ||||||
| Intercept | 3 | default | ||||||
| Intercept | 4 | default | ||||||
| Intercept | 5 | default | ||||||
| Intercept | 6 | default | ||||||
| Intercept | 7 | default | ||||||
| Intercept | 8 | default |
dataset1.brm.olr.priors = c(prior(normal(0, 1), class = "b"),
prior(student_t(3, 0, 1), class = "Intercept"))Step 1: Build the model (prior predictive checks)
TODO We need to fix this..
dataset1.brm.olr.priorchecks <- brm(
effectiveness ~ condition,
prior = dataset1.brm.olr.priors,
data = exp1.data,
family= cumulative("logit"),
backend = 'cmdstanr',
sample_prior = 'only',
file = "02-bayesian_stats/rds/dataset1.brm.bayesiant.priorchecks.rds"
)
dataset1.olr.yprior <-
posterior_predict(dataset1.brm.olr.priorchecks)
ggplot() +
geom_histogram(aes(x = dataset1.olr.yprior),
fill = '#351c75',
alpha = .5,
size = 1,
binwidth = .5) +
scale_x_continuous(breaks = 1:9, limits = c(1,9)) +
xlab('prior draws') + ylab('')
Step 2: Computing the posterior probability
dataset1.brm.olr1 =
brm(
dataset1.brm.olr.formula,
prior = dataset1.brm.olr.priors,
data = exp1.data,
backend = "cmdstanr",
file = "02-bayesian_stats/rds/dataset1.brm.olr1.rds"
)
summary(dataset1.brm.olr1)## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: effectiveness ~ condition
## Data: exp1.data (Number of observations: 123)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -4.56 0.96 -6.78 -2.99 1.00 1070 513
## Intercept[2] -3.99 0.78 -5.82 -2.71 1.00 1583 1198
## Intercept[3] -3.29 0.56 -4.55 -2.32 1.00 2496 2468
## Intercept[4] -2.14 0.34 -2.85 -1.50 1.00 3348 3016
## Intercept[5] -1.26 0.26 -1.79 -0.76 1.00 3758 3330
## Intercept[6] -0.27 0.23 -0.74 0.18 1.00 3652 3366
## Intercept[7] 1.26 0.26 0.77 1.81 1.00 3532 3334
## Intercept[8] 2.22 0.33 1.62 2.91 1.00 3937 2992
## conditionno_graph 0.18 0.31 -0.40 0.81 1.00 3434 2673
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 1.00 4000 4000
##
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).We see the error messages. Here we adjust a few sampling parameters to help the model converge.
dataset1.brm.olr2 =
brm(
dataset1.brm.olr.formula,
prior = dataset1.brm.olr.priors,
data = exp1.data,
backend = "cmdstanr",
warmup = 1500,
iter = 2500,
control = list(adapt_delta = 0.99, max_treedepth = 15),
file = "02-bayesian_stats/rds/dataset1.brm.olr2.rds"
)
summary(dataset1.brm.olr2)## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: effectiveness ~ condition
## Data: exp1.data (Number of observations: 123)
## Samples: 4 chains, each with iter = 2500; warmup = 1500; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -4.58 1.00 -6.88 -3.00 1.00 2194 1948
## Intercept[2] -3.99 0.78 -5.76 -2.71 1.00 2806 2384
## Intercept[3] -3.30 0.55 -4.48 -2.32 1.00 3705 2806
## Intercept[4] -2.15 0.34 -2.85 -1.51 1.00 4378 3270
## Intercept[5] -1.26 0.26 -1.78 -0.76 1.00 4525 3093
## Intercept[6] -0.27 0.24 -0.73 0.19 1.00 3935 3132
## Intercept[7] 1.27 0.26 0.76 1.79 1.00 4089 3216
## Intercept[8] 2.23 0.34 1.58 2.92 1.00 4890 3122
## conditionno_graph 0.19 0.31 -0.43 0.82 1.00 4856 2853
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 1.00 4000 4000
##
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Step 3: Evaluate the fit (MCMC traces)
plot(dataset1.brm.olr2)
Step 3: Evaluate posterior predictions
dataset1.olr.y <- exp1.data$effectiveness
dataset1.olr.yrep <- posterior_predict(dataset1.brm.olr2)
ppc_hist(y = dataset1.olr.y,
yrep = dataset1.olr.yrep[1000:1007, ],
binwidth = .5)
ppc_dens_overlay(y = dataset1.olr.y,
yrep = dataset1.olr.yrep[2000:2030, ])
ppc_stat_grouped(y = dataset1.olr.y,
yrep = dataset1.olr.yrep,
group = exp1.data$condition)
### Step 3: Compare posterior predictions
dataset1.olr.posterior_epred <-
tibble(condition = c('graph', 'no_graph')) %>%
add_epred_draws(dataset1.brm.olr2,
re_formula = NA,
allow_new_levels = FALSE) %>%
ungroup()
head(dataset1.olr.posterior_epred)| condition | .row | .chain | .iteration | .draw | .category | .epred |
|---|---|---|---|---|---|---|
| graph | 1 | NA | NA | 1 | 1 | 0.0045580 |
| graph | 1 | NA | NA | 2 | 1 | 0.0063148 |
| graph | 1 | NA | NA | 3 | 1 | 0.0069149 |
| graph | 1 | NA | NA | 4 | 1 | 0.0391315 |
| graph | 1 | NA | NA | 5 | 1 | 0.0081217 |
| graph | 1 | NA | NA | 6 | 1 | 0.0248435 |
dataset1.olr.posterior_comparison <-
dataset1.olr.posterior_epred %>%
select(-c(.chain, .iteration, .row)) %>%
group_by(.category) %>%
compare_levels(variable = .epred, by = condition)
head(dataset1.olr.posterior_comparison %>%
mean_qi())| .category | condition | .epred | .lower | .upper | .width | .point | .interval |
|---|---|---|---|---|---|---|---|
| 1 | no_graph - graph | -0.0025112 | -0.0174983 | 0.0069861 | 0.95 | mean | qi |
| 2 | no_graph - graph | -0.0013803 | -0.0110381 | 0.0040049 | 0.95 | mean | qi |
| 3 | no_graph - graph | -0.0026421 | -0.0162269 | 0.0074813 | 0.95 | mean | qi |
| 4 | no_graph - graph | -0.0102726 | -0.0506452 | 0.0227708 | 0.95 | mean | qi |
| 5 | no_graph - graph | -0.0137966 | -0.0621373 | 0.0329472 | 0.95 | mean | qi |
| 6 | no_graph - graph | -0.0145299 | -0.0680404 | 0.0330627 | 0.95 | mean | qi |
dataset1.olr.posterior_comparison %>%
mean_qi(.epred) %>%
ggplot() +
geom_point(aes(y = .epred, x = .category), size = 3, color = DARK_PURPLE) +
geom_errorbar(
aes(ymin = .lower, ymax = .upper, x = .category),
width = 0,
size = 2,
color = DARK_PURPLE,
alpha = .5
) +
geom_hline(aes(yintercept = 0), linetype = "dashed", color = "gray") +
xlab('') + ylab('no_graph - graph') +
theme(axis.title.y = element_text(angle = 0, vjust = .5))
Dataset 2
Load the data.
dataset2 = readr::read_csv("02-bayesian_stats/data/exp2.csv") %>%
mutate(condition = condition == 'expansive') %>%
group_by(participant)
head(dataset2)| participant | adjP10 | adjP20 | adjP30 | adjPDiff | condition | gender | index | change | subgroup | orig | BIS |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 56.50000 | 55.80000 | 58.40000 | 1.90000 | TRUE | female | 69 | 3.362832 | expansiveHigh | 56.90000 | high |
| 2 | 27.87500 | 40.14286 | 36.00000 | 8.12500 | FALSE | female | 60 | 29.147982 | constrictiveLow | 34.67262 | low |
| 3 | 61.00000 | 59.00000 | 76.50000 | 15.50000 | TRUE | female | 55 | 25.409836 | expansiveLow | 65.50000 | low |
| 4 | 24.57143 | 32.75000 | 37.85714 | 13.28571 | FALSE | female | 79 | 54.069767 | constrictiveHigh | 31.72619 | high |
| 5 | 64.71429 | 54.00000 | 41.00000 | -23.71429 | TRUE | female | 69 | -36.644592 | expansiveHigh | 53.23810 | high |
| 6 | 35.14286 | 52.40000 | 45.60000 | 10.45714 | FALSE | female | 43 | 29.756098 | constrictiveLow | 44.38095 | low |
We plot the data to give us a better picture about its distribution.
dataset2 %>%
mutate(c = as.factor(condition)) %>%
ggplot(aes(x = change)) +
geom_histogram(
aes(y = ..density..),
binwidth = 6,
fill = DARK_PURPLE,
color = 'white'
) +
geom_density(size = 1,
adjust = 1.5,
color = 'lightgray') +
scale_x_continuous(limits = c(-200, 200)) 
Step 1: Build the model (likelihood and priors)
This model is the BEST test model as described by Kruschke in the paper Bayesian estimation supersedes the t-test. In this model, \(\beta\) indicates the mean difference in the outcome variable between the two groups (in this case, the percent change in the BART scores). We fit different priors on \(\beta\) and set different weights on these priors to obtain our posterior estimate.
\[ \begin{align} y_{i} &\sim \mathrm{T}(\nu, \mu, \sigma) \\ \mu &= \alpha_{0} + \beta * x_i \\ \sigma &= \sigma_{a} + \sigma_{b}*x_i \\ \beta &\sim \mathrm{N}(\mu_{0}, \sigma_{0}) \\ \sigma_a, \sigma_b &\sim \mathrm{Cauchy}(0, 2) \\ \nu &\sim \mathrm{exp}(1/30) \end{align} \]
We translate the above specification to R code.
dataset2.brm.student.formula <- bf(change ~ condition,
sigma ~ condition,
family = student())
head(as_tibble(get_prior(dataset2.brm.student.formula, data = dataset2)))| prior | class | coef | group | resp | dpar | nlpar | bound | source |
|---|---|---|---|---|---|---|---|---|
| b | default | |||||||
| b | conditionTRUE | default | ||||||
| student_t(3, 22.6, 38.8) | Intercept | default | ||||||
| gamma(2, 0.1) | nu | default | ||||||
| b | sigma | default | ||||||
| b | conditionTRUE | sigma | default |
dataset2.brm.student.priorchecks = brm(
dataset2.brm.student.formula,
prior = c(
prior(normal(0, 2), class = "b"),
prior(cauchy(0, 2), class = "b", dpar = "sigma"),
prior(exponential(0.033), class = "nu"),
prior(student_t(3, 0, 5), class = "Intercept"),
prior(student_t(3, 0, 2), class = "Intercept", dpar = "sigma")
),
data = dataset2,
backend = "cmdstanr",
sample_prior = 'only',
file = "02-bayesian_stats/rds/dataset2.brm.student.priorchecks.rds"
)
# This breaks :-)
# dataset2.student.yprior <-
# posterior_predict(dataset2.brm.student.priorchecks)
# ggplot() +
# geom_density(
# aes(x = dataset2.student.yprior),
# color = '#351c75',
# alpha = .5,
# size = 1
# ) +
# xlab('prior draws') +
# ggtitle('prior preditive check')Step 2: Computing the posterior probability
dataset2.brm.student = brm(
bf(change ~ condition,
sigma ~ condition,
family = student()),
prior = c(
prior(normal(0, 2), class = "b"),
prior(cauchy(0, 2), class = "b", dpar = "sigma"),
prior(exponential(.033), class = "nu"),
prior(student_t(3, 0, 5), class = "Intercept"),
prior(student_t(3, 0, 2), class = "Intercept", dpar = "sigma")
),
data = dataset2,
backend = "cmdstanr",
file = "02-bayesian_stats/rds/dataset2.brm.student.rds"
)
dataset2.brm.student$prior| prior | class | coef | group | resp | dpar | nlpar | bound | source |
|---|---|---|---|---|---|---|---|---|
| normal(0, 2) | b | user | ||||||
| b | conditionTRUE | default | ||||||
| cauchy(0, 2) | b | sigma | user | |||||
| b | conditionTRUE | sigma | default | |||||
| student_t(3, 0, 5) | Intercept | user | ||||||
| student_t(3, 0, 2) | Intercept | sigma | user | |||||
| exponential(0.033) | nu | user |
Step 3: Evaluate the fit (summary)
summary(dataset2.brm.student)## Family: student
## Links: mu = identity; sigma = log; nu = identity
## Formula: change ~ condition
## sigma ~ condition
## Data: dataset2 (Number of observations: 80)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 21.15 4.95 11.41 31.04 1.00 2751 2071
## sigma_Intercept 3.37 0.20 2.98 3.76 1.00 2331 2293
## conditionTRUE -0.15 1.92 -4.01 3.46 1.00 3711 2961
## sigma_conditionTRUE 0.14 0.22 -0.31 0.55 1.00 3648 2440
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## nu 4.87 5.87 1.59 16.82 1.00 1781 1303
##
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Step 3: Evaluate the posterior predictions
plot(dataset2.brm.student)
dataset2.student.y <- dataset2.brm.student$data$change
dataset2.student.yrep <- posterior_predict(dataset2.brm.student)
dataset2.student.epred <- tibble(condition = c(TRUE, FALSE)) %>%
add_epred_draws(dataset2.brm.student,
re_formula = NA,
allow_new_levels = TRUE) %>%
ungroup()ppc_hist(y = dataset2.student.y,
yrep = dataset2.student.yrep[100:107, ],
binwidth = 10)
ppc_dens_overlay(y = dataset2.student.y,
yrep = dataset2.student.yrep[2000:2030, ])
ppc_stat_grouped(
y = dataset2.student.y,
yrep = dataset2.student.yrep,
group = dataset2$condition,
binwidth = 5
)
Step 3: Compare the posterior predictions
dataset2.student.epred_comparison <-
tibble(condition = c(TRUE, FALSE)) %>%
add_epred_draws(dataset2.brm.student,
re_formula = NA,
allow_new_levels = FALSE) %>%
ungroup() %>%
select(-c(.chain, .iteration, .row)) %>%
compare_levels(variable = .epred, by = condition) %>%
rename(diff = .epred)
kable(head(dataset2.student.epred_comparison))| .draw | condition | diff |
|---|---|---|
| 1 | TRUE - FALSE | -2.918310 |
| 2 | TRUE - FALSE | 2.367050 |
| 3 | TRUE - FALSE | 3.045840 |
| 4 | TRUE - FALSE | -0.455904 |
| 5 | TRUE - FALSE | -1.398820 |
| 6 | TRUE - FALSE | 2.193080 |
dataset2.student.epred_comparison %>%
select(diff) %>%
mean_qi() %>%
ggplot() +
geom_point(aes(x = diff, y = condition), size = 3, color = DARK_PURPLE) +
geom_errorbarh(
aes(xmin = .lower, xmax = .upper, y = condition),
height = 0,
color = DARK_PURPLE,
alpha = .5,
size = 2
) +
geom_vline(aes(xintercept = 0), linetype = "dashed", color = "gray") +
coord_cartesian(ylim = c(0, 2), xlim = c(-5, 5)) +
scale_y_discrete(label = c('expansive - not expansive')) +
xlab('') + ylab('')